03 - Clocks and cell counts

In this script, a trimmed-down version (excluding any plots that are the same as in Script 02, such as the loadings of the cell count PCs) of Script 02 will be run for several clocks. In other words, we will test how cell counts associate with several epigenetic clocks, both DNAmAge (the clock’s age estimation) and AgeAccel (the deviation of DNAmAge from calendar age). Since we see that the PC-approach works best, we will only apply this method for the clocks.

---

Setup

Load sample sheet.

setwd("/exports/molepi/RSC_BIOS/Users/tjonkman/cellcounts")
.libPaths("/exports/molepi/RSC_BIOS/Users/tjonkman/Packages/4.3.1")

load("01_sample_sheet.rda")
dim(ss)

## [1] 4058   27

#Load PCA data.
load("02_pca_df.rda")
ss.pc <- cbind(ss, pca.df)

library(ggplot2)
library(reshape2)

#Print the correlations between the investigated clocks (residuals for all except dunedinPACE).
round(cor(ss[,c(5:10)]), 2)

##             hannum horvath zhang phenoage grimage dunedinpace
## hannum        1.00    0.96  0.98     0.96    0.95        0.49
## horvath       0.96    1.00  0.97     0.95    0.93        0.45
## zhang         0.98    0.97  1.00     0.96    0.96        0.45
## phenoage      0.96    0.95  0.96     1.00    0.95        0.54
## grimage       0.95    0.93  0.96     0.95    1.00        0.57
## dunedinpace   0.49    0.45  0.45     0.54    0.57        1.00

round(cor(ss[,c(23:27, 10)]), 2)

##             hannumRes horvathRes zhangRes phenoageRes grimageRes dunedinpace
## hannumRes        1.00       0.55     0.66        0.59       0.39        0.27
## horvathRes       0.55       1.00     0.60        0.55       0.22        0.15
## zhangRes         0.66       0.60     1.00        0.55       0.23        0.15
## phenoageRes      0.59       0.55     0.55        1.00       0.49        0.40
## grimageRes       0.39       0.22     0.23        0.49       1.00        0.60
## dunedinpace      0.27       0.15     0.15        0.40       0.60        1.00

---

Correlation between DNAmAge and calendar age

#DNAmAge
plot.data <- ss[,c("age", "hannum", "horvath", "zhang", "phenoage", "grimage")]
colnames(plot.data) <- c("Calendar Age", "Hannum", "Horvath", "Zhang", "PhenoAge", "GrimAge")
plot.data <- melt(plot.data, id.vars = "Calendar Age", variable.name = "Clock", value.name = "DNAmAge")

library(ggpubr)
ggplot(data = plot.data, aes(x = `Calendar Age`, y = DNAmAge)) + 
  geom_point(shape = 1, size = 1) + 
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", color= "red") +
  geom_smooth(method = "lm") +
  stat_cor(vjust = 0.5) +
  facet_wrap(facets = vars(Clock), ncol = 2, scales = "free", dir = "v") +
  scale_x_continuous(limits = c(18, 87), breaks = c(0, 20, 40, 60, 80)) +
  scale_y_continuous(limits = c(-5, 123), breaks = c(0, 20, 40, 60, 80, 100, 120)) +
  theme_bw()

## `geom_smooth()` using formula = 'y ~ x'

#AgeAccel
plot.data <- ss[,c("age", "hannumRes", "horvathRes", "zhangRes", "phenoageRes", "grimageRes")]
colnames(plot.data) <- c("Calendar Age", "Hannum", "Horvath", "Zhang", "PhenoAge", "GrimAge")
plot.data <- melt(plot.data, id.vars = "Calendar Age", variable.name = "Clock", value.name = "AgeAccel")

ggplot(data = plot.data, aes(x = `Calendar Age`, y = AgeAccel)) + 
  geom_point(shape = 1, size = 1) + 
  geom_abline(slope = 0, intercept = 0, linetype = "dashed", color= "red") +
  geom_smooth(method = "lm") +
  stat_cor(vjust = 0.5) +
  facet_wrap(facets = vars(Clock), ncol = 2, scales = "free", dir = "v") + 
  scale_x_continuous(limits = c(18, 87), breaks = c(0, 20, 40, 60, 80)) +
  scale_y_continuous(limits = c(-38, 68), breaks = c(-20, 0, 20, 40, 60)) +
  theme_bw()

## `geom_smooth()` using formula = 'y ~ x'

#DunedinPACE
plot.data <- ss[,c("age", "dunedinpace")]
colnames(plot.data) <- c("Calendar Age", "DunedinPACE")
plot.data <- melt(plot.data, id.vars = "Calendar Age", variable.name = "Clock", value.name = "AgeAccel")

ggplot(data = plot.data, aes(x = `Calendar Age`, y = AgeAccel)) + 
  geom_point(shape = 1, size = 1) + 
  geom_abline(slope = 0, intercept = 52, linetype = "dashed", color = "red") +
  geom_smooth(method = "lm", ) +
  stat_cor(vjust = 0.5) +
  facet_wrap(facets = vars(Clock), nrow = 1, scales = "free") + 
  scale_x_continuous(limits = c(18, 87), breaks = c(0, 20, 40, 60, 80)) +
  scale_y_continuous(limits = c(26, 86), breaks = c(20, 40, 60, 80)) +
  theme_bw()

## `geom_smooth()` using formula = 'y ~ x'

---

Clock ~ PC models

Run a clock ~ PCs model for each clock’s DNAmAge and AgeAccel.

var.df <- data.frame(
  row.names = c("age", "ageRes", "hannum", "hannumRes", "horvath", "horvathRes", "zhang", "zhangRes", "phenoage", "phenoageRes", "grimage", "grimageRes", "dunedinpaceAge", "dunedinpace"),
  variable = c("Age", "AgeAccel", "Hannum", "Hannum residual", "Horvath", "Horvath residual", "Zhang", "Zhang residual", "PhenoAge", "PhenoAge residual", "GrimAge", "GrimAge residual", "DunedinPACE-age", "DunedinPACE"),
  var.ex = NaN
)

#Prepare PC labels.
pc.labs <- c("PC1: Neutrophils", "PC2: T cell naive/memory", "PC3: Treg/CD4Tmem", "PC4: Bmem/Mono", "PC5: Bnv/Eos", "PC6: Naive CD8 T cells", "PC7: Mixed (CD4Tnv/CD8Tnv)", "PC8: Eos/Mono", "PC9: Mixed (CD4Tnv+CD8Tmem)", "PC10: Mixed (Bmem)", "PC11: Mixed (Baso)")

PCEffects <- function(form){

    fit <- lm(formula = form, data = ss.pc)
    var.ex <- round(summary(fit)$r.squared * 100, 1)
    print(paste0("Variance explained: ", var.ex, "%"))
    
    #Store variance explained.
    var.df[as.character(form)[2], "var.ex"] <<- var.ex
    
    res <- as.data.frame(summary(fit)$coefficients)
    colnames(res) <- c("estimate", "std.error", "t.stat", "p.val")
    
    #Calculate 95% confidence intervals.
    res$ci.lower <- res$estimate - (1.96 * res$std.error)
    res$ci.upper <- res$estimate + (1.96 * res$std.error)
    
    #Select only cell type effects, removing the intercept.
    res <- res[(rownames(res) %in% c("PC1", "PC2", "PC3", "PC4", "PC5", "PC6", "PC7", "PC8", "PC9", "PC10", "PC11")),]
    
    res$PC <- factor(rownames(res), levels = rownames(res))
    
    #Add PC labels.
    res$PC.label <- pc.labs
    res$PC.label <- factor(res$PC.label, levels = unique(res$PC.label))
    
    print(res)
    return(res)
}

res.age <- PCEffects(form = formula(age ~ PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC7 + PC8 + PC9 + PC10 + PC11))

## [1] "Variance explained: 43.5%"
##         estimate  std.error      t.stat         p.val   ci.lower    ci.upper
## PC1   0.02332937 0.09791573   0.2382597  8.116918e-01 -0.1685855  0.21524419
## PC2  -6.46156949 0.14944642 -43.2366974  0.000000e+00 -6.7544845 -6.16865452
## PC3  -0.72824707 0.17574455  -4.1437817  3.486108e-05 -1.0727064 -0.38378774
## PC4   0.99617076 0.19500859   5.1083429  3.399769e-07  0.6139539  1.37838760
## PC5  -1.95336043 0.22412537  -8.7154809  4.155524e-18 -2.3926462 -1.51407469
## PC6  -6.28966305 0.22680381 -27.7317344 4.057491e-155 -6.7341985 -5.84512758
## PC7   4.35027432 0.23649668  18.3946525  1.229801e-72  3.8867408  4.81380783
## PC8  -0.44156734 0.24144279  -1.8288695  6.749274e-02 -0.9147952  0.03166052
## PC9   0.83517919 0.26369586   3.1672063  1.550561e-03  0.3183353  1.35202307
## PC10  0.09467097 0.29267373   0.3234693  7.463566e-01 -0.4789695  0.66831147
## PC11  0.76572771 0.31814402   2.4068587  1.613492e-02  0.1421654  1.38928999
##        PC                    PC.label
## PC1   PC1            PC1: Neutrophils
## PC2   PC2    PC2: T cell naive/memory
## PC3   PC3           PC3: Treg/CD4Tmem
## PC4   PC4              PC4: Bmem/Mono
## PC5   PC5                PC5: Bnv/Eos
## PC6   PC6      PC6: Naive CD8 T cells
## PC7   PC7  PC7: Mixed (CD4Tnv/CD8Tnv)
## PC8   PC8               PC8: Eos/Mono
## PC9   PC9 PC9: Mixed (CD4Tnv+CD8Tmem)
## PC10 PC10          PC10: Mixed (Bmem)
## PC11 PC11          PC11: Mixed (Baso)

#Clocks.
res.hannum <- PCEffects(form = formula(hannum ~ PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC7 + PC8 + PC9 + PC10 + PC11))

## [1] "Variance explained: 52.6%"
##        estimate  std.error     t.stat         p.val   ci.lower    ci.upper   PC
## PC1   0.6132778 0.08875804   6.909547  5.623812e-12  0.4393120  0.78724355  PC1
## PC2  -7.2567293 0.13546926 -53.567351  0.000000e+00 -7.5222490 -6.99120955  PC2
## PC3  -0.3002849 0.15930783  -1.884935  5.951021e-02 -0.6125283  0.01195843  PC3
## PC4   0.8487059 0.17677018   4.801183  1.634319e-06  0.5022364  1.19517544  PC4
## PC5  -2.0269365 0.20316377  -9.976860  3.557852e-23 -2.4251375 -1.62873548  PC5
## PC6  -6.8969005 0.20559171 -33.546589 6.572806e-218 -7.2998602 -6.49394076  PC6
## PC7   3.6262222 0.21437804  16.915082  4.521235e-62  3.2060412  4.04640313  PC7
## PC8  -0.4476190 0.21886156  -2.045215  4.089820e-02 -0.8765876 -0.01865032  PC8
## PC9   0.6978350 0.23903338   2.919404  3.526367e-03  0.2293296  1.16634047  PC9
## PC10  0.3819960 0.26530106   1.439859  1.499847e-01 -0.1379940  0.90198611 PC10
## PC11  1.2517551 0.28838921   4.340506  1.456290e-05  0.6865122  1.81699795 PC11
##                         PC.label
## PC1             PC1: Neutrophils
## PC2     PC2: T cell naive/memory
## PC3            PC3: Treg/CD4Tmem
## PC4               PC4: Bmem/Mono
## PC5                 PC5: Bnv/Eos
## PC6       PC6: Naive CD8 T cells
## PC7   PC7: Mixed (CD4Tnv/CD8Tnv)
## PC8                PC8: Eos/Mono
## PC9  PC9: Mixed (CD4Tnv+CD8Tmem)
## PC10          PC10: Mixed (Bmem)
## PC11          PC11: Mixed (Baso)

res.horvath <- PCEffects(form = formula(horvath ~ PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC7 + PC8 + PC9 + PC10 + PC11))

## [1] "Variance explained: 43.8%"
##        estimate  std.error     t.stat         p.val    ci.lower   ci.upper   PC
## PC1   0.1219737 0.09108413   1.339132  1.806029e-01 -0.05655122  0.3004986  PC1
## PC2  -6.4140141 0.13901952 -46.137507  0.000000e+00 -6.68649232 -6.1415358  PC2
## PC3  -0.4607825 0.16348283  -2.818538  4.847731e-03 -0.78120888 -0.1403562  PC3
## PC4   0.7518009 0.18140282   4.144373  3.477167e-05  0.39625138  1.1073504  PC4
## PC5  -1.5033864 0.20848812  -7.210897  6.601465e-13 -1.91202308 -1.0947497  PC5
## PC6  -5.4686838 0.21097968 -25.920429 3.519637e-137 -5.88220393 -5.0551636  PC6
## PC7   3.3604922 0.21999628  15.275223  2.937058e-51  2.92929951  3.7916849  PC7
## PC8  -0.5768126 0.22459729  -2.568208  1.025815e-02 -1.01702329 -0.1366019  PC8
## PC9   1.0675295 0.24529777   4.351974  1.382481e-05  0.58674585  1.5483131  PC9
## PC10  0.1074964 0.27225384   0.394839  6.929825e-01 -0.42612111  0.6411140 PC10
## PC11  1.0633693 0.29594707   3.593106  3.306486e-04  0.48331303  1.6434256 PC11
##                         PC.label
## PC1             PC1: Neutrophils
## PC2     PC2: T cell naive/memory
## PC3            PC3: Treg/CD4Tmem
## PC4               PC4: Bmem/Mono
## PC5                 PC5: Bnv/Eos
## PC6       PC6: Naive CD8 T cells
## PC7   PC7: Mixed (CD4Tnv/CD8Tnv)
## PC8                PC8: Eos/Mono
## PC9  PC9: Mixed (CD4Tnv+CD8Tmem)
## PC10          PC10: Mixed (Bmem)
## PC11          PC11: Mixed (Baso)

res.zhang <- PCEffects(form = formula(zhang ~ PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC7 + PC8 + PC9 + PC10 + PC11))

## [1] "Variance explained: 46%"
##         estimate  std.error      t.stat         p.val    ci.lower    ci.upper
## PC1   0.10182666 0.09331251   1.0912434  2.752308e-01 -0.08106585  0.28471918
## PC2  -6.54800091 0.14242064 -45.9764891  0.000000e+00 -6.82714536 -6.26885646
## PC3  -0.62004035 0.16748245  -3.7021214  2.166524e-04 -0.94830594 -0.29177475
## PC4   0.90703455 0.18584084   4.8807061  1.097805e-06  0.54278649  1.27128260
## PC5  -1.95795466 0.21358879  -9.1669356  7.560872e-20 -2.37658869 -1.53932064
## PC6  -6.35852571 0.21614130 -29.4183739 1.488056e-172 -6.78216267 -5.93488876
## PC7   4.05406403 0.22537850  17.9878032  1.174807e-69  3.61232218  4.49580588
## PC8  -0.48414463 0.23009208  -2.1041343  3.542828e-02 -0.93512509 -0.03316416
## PC9   0.90632251 0.25129898   3.6065506  3.140400e-04  0.41377650  1.39886851
## PC10  0.07382008 0.27891454   0.2646692  7.912778e-01 -0.47285242  0.62049258
## PC11  0.74454187 0.30318743   2.4557149  1.410223e-02  0.15029451  1.33878923
##        PC                    PC.label
## PC1   PC1            PC1: Neutrophils
## PC2   PC2    PC2: T cell naive/memory
## PC3   PC3           PC3: Treg/CD4Tmem
## PC4   PC4              PC4: Bmem/Mono
## PC5   PC5                PC5: Bnv/Eos
## PC6   PC6      PC6: Naive CD8 T cells
## PC7   PC7  PC7: Mixed (CD4Tnv/CD8Tnv)
## PC8   PC8               PC8: Eos/Mono
## PC9   PC9 PC9: Mixed (CD4Tnv+CD8Tmem)
## PC10 PC10          PC10: Mixed (Bmem)
## PC11 PC11          PC11: Mixed (Baso)

res.phenoage <- PCEffects(form = formula(phenoage ~ PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC7 + PC8 + PC9 + PC10 + PC11))

## [1] "Variance explained: 49.7%"
##         estimate  std.error     t.stat         p.val    ci.lower    ci.upper
## PC1   1.00896092 0.09640441  10.465920  2.600157e-25  0.82000827  1.19791357
## PC2  -7.25540535 0.14713974 -49.309626  0.000000e+00 -7.54379924 -6.96701147
## PC3  -0.24299304 0.17303197  -1.404325  1.602990e-01 -0.58213569  0.09614962
## PC4   0.94273514 0.19199867   4.910113  9.461544e-07  0.56641775  1.31905253
## PC5  -2.53284680 0.22066604 -11.478190  4.929291e-30 -2.96535224 -2.10034137
## PC6  -7.16087308 0.22330313 -32.067947 2.970320e-201 -7.59854723 -6.72319894
## PC7   3.73178703 0.23284640  16.026819  4.280387e-56  3.27540809  4.18816598
## PC8  -0.76464931 0.23771616  -3.216648  1.307192e-03 -1.23057299 -0.29872563
## PC9   0.27279119 0.25962576   1.050709  2.934549e-01 -0.23607530  0.78165768
## PC10  0.08369212 0.28815636   0.290440  7.714946e-01 -0.48109434  0.64847859
## PC11  0.54496334 0.31323353   1.739799  8.197037e-02 -0.06897437  1.15890105
##        PC                    PC.label
## PC1   PC1            PC1: Neutrophils
## PC2   PC2    PC2: T cell naive/memory
## PC3   PC3           PC3: Treg/CD4Tmem
## PC4   PC4              PC4: Bmem/Mono
## PC5   PC5                PC5: Bnv/Eos
## PC6   PC6      PC6: Naive CD8 T cells
## PC7   PC7  PC7: Mixed (CD4Tnv/CD8Tnv)
## PC8   PC8               PC8: Eos/Mono
## PC9   PC9 PC9: Mixed (CD4Tnv+CD8Tmem)
## PC10 PC10          PC10: Mixed (Bmem)
## PC11 PC11          PC11: Mixed (Baso)

res.grimage <- PCEffects(form = formula(grimage ~ PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC7 + PC8 + PC9 + PC10 + PC11))

## [1] "Variance explained: 46.6%"
##         estimate  std.error      t.stat         p.val    ci.lower   ci.upper
## PC1   0.72059896 0.07978605   9.0316407  2.562624e-19  0.56421829  0.8769796
## PC2  -5.57424622 0.12177553 -45.7747638  0.000000e+00 -5.81292626 -5.3355662
## PC3  -0.59797093 0.14320442  -4.1756459  3.033989e-05 -0.87865158 -0.3172903
## PC4   1.08576990 0.15890160   6.8329701  9.561125e-12  0.77432276  1.3972170
## PC5  -2.36538064 0.18262724 -12.9519595  1.269433e-37 -2.72333003 -2.0074312
## PC6  -5.58907745 0.18480975 -30.2423301 2.588131e-181 -5.95130456 -5.2268503
## PC7   2.64687473 0.19270793  13.7351624  5.398072e-42  2.26916718  3.0245823
## PC8  -0.91678279 0.19673823  -4.6599116  3.265219e-06 -1.30238973 -0.5311759
## PC9   0.02666734 0.21487102   0.1241086  9.012355e-01 -0.39447986  0.4478145
## PC10 -0.12448889 0.23848346  -0.5220022  6.016974e-01 -0.59191648  0.3429387
## PC11  0.47650861 0.25923778   1.8381141  6.611887e-02 -0.03159744  0.9846147
##        PC                    PC.label
## PC1   PC1            PC1: Neutrophils
## PC2   PC2    PC2: T cell naive/memory
## PC3   PC3           PC3: Treg/CD4Tmem
## PC4   PC4              PC4: Bmem/Mono
## PC5   PC5                PC5: Bnv/Eos
## PC6   PC6      PC6: Naive CD8 T cells
## PC7   PC7  PC7: Mixed (CD4Tnv/CD8Tnv)
## PC8   PC8               PC8: Eos/Mono
## PC9   PC9 PC9: Mixed (CD4Tnv+CD8Tmem)
## PC10 PC10          PC10: Mixed (Bmem)
## PC11 PC11          PC11: Mixed (Baso)

#Clock residuals (and also DunedinPACE).
res.hannumRes <- PCEffects(form = formula(hannumRes ~ PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC7 + PC8 + PC9 + PC10 + PC11))

## [1] "Variance explained: 21.2%"
##         estimate  std.error      t.stat         p.val    ci.lower    ci.upper
## PC1   0.59118817 0.03324889  17.7806917  3.676070e-68  0.52602036  0.65635599
## PC2  -1.13853297 0.05074697 -22.4354852 3.607497e-105 -1.23799704 -1.03906890
## PC3   0.38926249 0.05967694   6.5228297  7.749891e-11  0.27229570  0.50622929
## PC4  -0.09452752 0.06621836  -1.4275123  1.535094e-01 -0.22431550  0.03526046
## PC5  -0.17737925 0.07610543  -2.3307042  1.981786e-02 -0.32654590 -0.02821260
## PC6  -0.94147537 0.07701494 -12.2245808  8.996195e-34 -1.09242466 -0.79052609
## PC7  -0.49287496 0.08030631  -6.1374371  9.191863e-10 -0.65027534 -0.33547458
## PC8  -0.02951689 0.08198585  -0.3600242  7.188478e-01 -0.19020915  0.13117537
## PC9  -0.09296203 0.08954224  -1.0381919  2.992428e-01 -0.26846481  0.08254076
## PC10  0.29235596 0.09938215   2.9417352  3.282309e-03  0.09756695  0.48714497
## PC11  0.52671879 0.10803100   4.8756262  1.126269e-06  0.31497802  0.73845956
##        PC                    PC.label
## PC1   PC1            PC1: Neutrophils
## PC2   PC2    PC2: T cell naive/memory
## PC3   PC3           PC3: Treg/CD4Tmem
## PC4   PC4              PC4: Bmem/Mono
## PC5   PC5                PC5: Bnv/Eos
## PC6   PC6      PC6: Naive CD8 T cells
## PC7   PC7  PC7: Mixed (CD4Tnv/CD8Tnv)
## PC8   PC8               PC8: Eos/Mono
## PC9   PC9 PC9: Mixed (CD4Tnv+CD8Tmem)
## PC10 PC10          PC10: Mixed (Bmem)
## PC11 PC11          PC11: Mixed (Baso)

res.horvathRes <- PCEffects(form = formula(horvathRes ~ PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC7 + PC8 + PC9 + PC10 + PC11))

## [1] "Variance explained: 5.4%"
##         estimate  std.error      t.stat        p.val    ci.lower     ci.upper
## PC1   0.10132432 0.03720756   2.7232186 6.492667e-03  0.02839750  0.174251133
## PC2  -0.69473001 0.05678900 -12.2335307 8.089432e-34 -0.80603645 -0.583423561
## PC3   0.18380574 0.06678218   2.7523170 5.943903e-03  0.05291266  0.314698818
## PC4  -0.12993279 0.07410244  -1.7534213 7.960543e-02 -0.27517356  0.015307984
## PC5   0.22557796 0.08516669   2.6486642 8.112449e-03  0.05865125  0.392504670
## PC6   0.09844194 0.08618448   1.1422235 2.534287e-01 -0.07047964  0.267363531
## PC7  -0.49003583 0.08986773  -5.4528562 5.252568e-08 -0.66617659 -0.313895071
## PC8  -0.18597116 0.09174723  -2.0269947 4.272840e-02 -0.36579574 -0.006146583
## PC9   0.32829313 0.10020331   3.2762704 1.060748e-03  0.13189465  0.524691603
## PC10  0.02370097 0.11121477   0.2131099 8.312520e-01 -0.19427998  0.241681923
## PC11  0.38560605 0.12089338   3.1896376 1.435391e-03  0.14865504  0.622557070
##        PC                    PC.label
## PC1   PC1            PC1: Neutrophils
## PC2   PC2    PC2: T cell naive/memory
## PC3   PC3           PC3: Treg/CD4Tmem
## PC4   PC4              PC4: Bmem/Mono
## PC5   PC5                PC5: Bnv/Eos
## PC6   PC6      PC6: Naive CD8 T cells
## PC7   PC7  PC7: Mixed (CD4Tnv/CD8Tnv)
## PC8   PC8               PC8: Eos/Mono
## PC9   PC9 PC9: Mixed (CD4Tnv+CD8Tmem)
## PC10 PC10          PC10: Mixed (Bmem)
## PC11 PC11          PC11: Mixed (Baso)

res.zhangRes <- PCEffects(form = formula(zhangRes ~ PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC7 + PC8 + PC9 + PC10 + PC11))

## [1] "Variance explained: 4%"
##         estimate  std.error      t.stat        p.val     ci.lower     ci.upper
## PC1   0.07948475 0.02331971   3.4084791 6.596427e-04  0.033778116  0.125191387
## PC2  -0.35992888 0.03559232 -10.1125445 9.291185e-24 -0.429689814 -0.290167937
## PC3   0.07738236 0.04185551   1.8487976 6.455995e-02 -0.004654432  0.159419158
## PC4  -0.04697150 0.04644345  -1.0113698 3.119000e-01 -0.138000668  0.044057661
## PC5  -0.08727370 0.05337793  -1.6350146 1.021238e-01 -0.191894454  0.017347049
## PC6  -0.33508387 0.05401583  -6.2034379 6.078539e-10 -0.440954906 -0.229212838
## PC7  -0.11207717 0.05632430  -1.9898547 4.667414e-02 -0.222472796 -0.001681545
## PC8  -0.06126741 0.05750227  -1.0654781 2.867232e-01 -0.173971855  0.051437036
## PC9   0.10649377 0.06280208   1.6957044 9.001881e-02 -0.016598315  0.229585850
## PC10 -0.01684377 0.06970348  -0.2416489 8.090645e-01 -0.153462596  0.119775054
## PC11  0.01122495 0.07576951   0.1481461 8.822349e-01 -0.137283294  0.159733203
##        PC                    PC.label
## PC1   PC1            PC1: Neutrophils
## PC2   PC2    PC2: T cell naive/memory
## PC3   PC3           PC3: Treg/CD4Tmem
## PC4   PC4              PC4: Bmem/Mono
## PC5   PC5                PC5: Bnv/Eos
## PC6   PC6      PC6: Naive CD8 T cells
## PC7   PC7  PC7: Mixed (CD4Tnv/CD8Tnv)
## PC8   PC8               PC8: Eos/Mono
## PC9   PC9 PC9: Mixed (CD4Tnv+CD8Tmem)
## PC10 PC10          PC10: Mixed (Bmem)
## PC11 PC11          PC11: Mixed (Baso)

res.phenoageRes <- PCEffects(form = formula(phenoageRes ~ PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC7 + PC8 + PC9 + PC10 + PC11))

## [1] "Variance explained: 19.9%"
##          estimate  std.error       t.stat         p.val   ci.lower    ci.upper
## PC1   0.986125092 0.04216812  23.38556170 1.492941e-113  0.9034756  1.06877460
## PC2  -0.930532010 0.06436018 -14.45819420  3.105188e-46 -1.0566780 -0.80438605
## PC3   0.469847779 0.07568567   6.20788328  5.910684e-10  0.3215039  0.61819168
## PC4  -0.032361369 0.08398186  -0.38533759  7.000076e-01 -0.1969658  0.13224308
## PC5  -0.620810215 0.09652122  -6.43185237  1.407249e-10 -0.8099918 -0.43162863
## PC6  -1.004269455 0.09767470 -10.28177625  1.700205e-24 -1.1957119 -0.81282703
## PC7  -0.526456128 0.10184901  -5.16898618  2.467332e-07 -0.7260802 -0.32683207
## PC8  -0.332423434 0.10397909  -3.19702206  1.399237e-03 -0.5362224 -0.12862442
## PC9  -0.544719574 0.11356253  -4.79664871  1.671518e-06 -0.7673021 -0.32213701
## PC10 -0.008976056 0.12604206  -0.07121476  9.432304e-01 -0.2560185  0.23806638
## PC11 -0.204565200 0.13701103  -1.49305645  1.355004e-01 -0.4731068  0.06397642
##        PC                    PC.label
## PC1   PC1            PC1: Neutrophils
## PC2   PC2    PC2: T cell naive/memory
## PC3   PC3           PC3: Treg/CD4Tmem
## PC4   PC4              PC4: Bmem/Mono
## PC5   PC5                PC5: Bnv/Eos
## PC6   PC6      PC6: Naive CD8 T cells
## PC7   PC7  PC7: Mixed (CD4Tnv/CD8Tnv)
## PC8   PC8               PC8: Eos/Mono
## PC9   PC9 PC9: Mixed (CD4Tnv+CD8Tmem)
## PC10 PC10          PC10: Mixed (Bmem)
## PC11 PC11          PC11: Mixed (Baso)

res.grimageRes <- PCEffects(form = formula(grimageRes ~ PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC7 + PC8 + PC9 + PC10 + PC11))

## [1] "Variance explained: 28.3%"
##          estimate  std.error      t.stat         p.val    ci.lower    ci.upper
## PC1   0.701748692 0.02471451  28.3942038 6.940667e-162  0.65330826  0.75018912
## PC2  -0.353260857 0.03772115  -9.3650595  1.227348e-20 -0.42719432 -0.27932739
## PC3  -0.009543107 0.04435896  -0.2151337  8.296739e-01 -0.09648667  0.07740045
## PC4   0.280858194 0.04922132   5.7060277  1.239179e-08  0.18438441  0.37733197
## PC5  -0.787054164 0.05657056 -13.9127865  5.112909e-43 -0.89793247 -0.67617586
## PC6  -0.506993480 0.05724662  -8.8563048  1.215313e-18 -0.61919685 -0.39479011
## PC7  -0.868171979 0.05969315 -14.5439120  9.473587e-47 -0.98517056 -0.75117340
## PC8  -0.559993834 0.06094158  -9.1890271  6.184599e-20 -0.67943933 -0.44054834
## PC9  -0.648162255 0.06655839  -9.7382508  3.620692e-22 -0.77861669 -0.51770782
## PC10 -0.200983577 0.07387257  -2.7206794  6.542653e-03 -0.34577381 -0.05619334
## PC11 -0.142203784 0.08030142  -1.7708750  7.665676e-02 -0.29959457  0.01518700
##        PC                    PC.label
## PC1   PC1            PC1: Neutrophils
## PC2   PC2    PC2: T cell naive/memory
## PC3   PC3           PC3: Treg/CD4Tmem
## PC4   PC4              PC4: Bmem/Mono
## PC5   PC5                PC5: Bnv/Eos
## PC6   PC6      PC6: Naive CD8 T cells
## PC7   PC7  PC7: Mixed (CD4Tnv/CD8Tnv)
## PC8   PC8               PC8: Eos/Mono
## PC9   PC9 PC9: Mixed (CD4Tnv+CD8Tmem)
## PC10 PC10          PC10: Mixed (Bmem)
## PC11 PC11          PC11: Mixed (Baso)

res.dunedinpace <- PCEffects(form = formula(dunedinpace ~ PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC7 + PC8 + PC9 + PC10 + PC11))

## [1] "Variance explained: 26.3%"
##         estimate  std.error       t.stat        p.val   ci.lower   ci.upper
## PC1   0.92892594 0.04476168  20.75270619 5.446016e-91  0.8411931  1.0166588
## PC2  -1.45073602 0.06831867 -21.23483952 5.796187e-95 -1.5846406 -1.3168314
## PC3  -0.28795801 0.08034073  -3.58420945 3.420878e-04 -0.4454258 -0.1304902
## PC4   0.57841205 0.08914719   6.48828144 9.729261e-11  0.4036836  0.7531405
## PC5  -0.48324924 0.10245778  -4.71656960 2.479421e-06 -0.6840665 -0.2824320
## PC6  -2.11675678 0.10368221 -20.41581422 2.942103e-88 -2.3199739 -1.9135396
## PC7  -0.08865388 0.10811326  -0.82000929 4.122591e-01 -0.3005559  0.1232481
## PC8  -0.86827802 0.11037435  -7.86666512 4.638406e-15 -1.0846117 -0.6519443
## PC9  -0.30841626 0.12054722  -2.55846835 1.054953e-02 -0.5446888 -0.0721437
## PC10  0.01309479 0.13379431   0.09787256 9.220383e-01 -0.2491421  0.2753316
## PC11 -0.14563818 0.14543793  -1.00137690 3.167045e-01 -0.4306965  0.1394202
##        PC                    PC.label
## PC1   PC1            PC1: Neutrophils
## PC2   PC2    PC2: T cell naive/memory
## PC3   PC3           PC3: Treg/CD4Tmem
## PC4   PC4              PC4: Bmem/Mono
## PC5   PC5                PC5: Bnv/Eos
## PC6   PC6      PC6: Naive CD8 T cells
## PC7   PC7  PC7: Mixed (CD4Tnv/CD8Tnv)
## PC8   PC8               PC8: Eos/Mono
## PC9   PC9 PC9: Mixed (CD4Tnv+CD8Tmem)
## PC10 PC10          PC10: Mixed (Bmem)
## PC11 PC11          PC11: Mixed (Baso)

---

Plot PC effect sizes

Plot the effect sizes of the clocks versus the effect sizes of the residuals.

plotDNAmAgeAndResidual <- function(dnamage, age.residual, clock.lab){
  
  dat <- rbind(dnamage, age.residual)
  dat$type <- factor(c(rep("DNAmAge", 11), rep("AgeAccel", 11)), levels = c("DNAmAge", "AgeAccel"))
  dat
  
  #Prepare annotations for variance explained for the clock and its residual.
  clock.var <- paste0("Variance explained: ", round(var.df[grep(clock.lab, var.df$variable)[1],"var.ex"], 1), "%")
  resid.var <- paste0("Variance explained: ", round(var.df[grep(clock.lab, var.df$variable)[2],"var.ex"], 1), "%")
  
  p <- ggplot(dat, aes(x = PC.label, y = estimate)) +
      geom_hline(yintercept = 0, linetype = "dashed") +
      geom_point(aes(color = type), shape = 19, size = 1.5) +
      geom_errorbar(aes(ymin = ci.lower, ymax = ci.upper, color = type), width = 0.4) +
      labs(x = NULL, y = paste0(clock.lab, " effect"), color = "") +
      theme_bw() +
      scale_color_manual(values = c("#005580", "#993300"), drop = F) +
      scale_y_continuous(limits = c(-7.6, 4.7), breaks = seq(-10, 10, 2)) +
      theme(axis.text.x = element_text(angle = 45, hjust = 1),
            plot.margin = margin(t = 5.5, r = 5.5, b = 5.5, l = 22, unit = "pt"))
  print(p)
  
}

plotDNAmAgeAndResidual(res.hannum, res.hannumRes, "Hannum")

plotDNAmAgeAndResidual(res.horvath, res.horvathRes, "Horvath")

plotDNAmAgeAndResidual(res.zhang, res.zhangRes, "Zhang")

plotDNAmAgeAndResidual(res.phenoage, res.phenoageRes, "PhenoAge")

plotDNAmAgeAndResidual(res.grimage, res.grimageRes, "GrimAge")

#For DunedinPACE, we only have AgeAccel (the clock itself is a putative age acceleration marker). Add a dummy dataframe for the DNAmAge.
res.dunedinpace.dummy <- res.dunedinpace
for(i in 1:7){
  if(class(res.dunedinpace.dummy[,i]) == "numeric"){
    res.dunedinpace.dummy[,i] <- NaN
  } else{
    res.dunedinpace.dummy[,i] <- NA
  }
}

plotDNAmAgeAndResidual(res.dunedinpace.dummy, res.dunedinpace, "DunedinPACE")

## Warning: Removed 11 rows containing missing values (`geom_point()`).

Make the plots for all clocks in a single faceted figure.

plot.data <- rbind(res.hannum, res.hannumRes, res.horvath, res.horvathRes, res.zhang, res.zhangRes, res.phenoage, res.phenoageRes, res.grimage, res.grimageRes, res.dunedinpace.dummy, res.dunedinpace)
plot.data$type <- factor(rep(c(rep("DNAmAge", 11), rep("AgeAccel", 11)), 6), levels = c("DNAmAge", "AgeAccel"))
plot.data$clock <- factor(c(rep("Hannum", 22), rep("Horvath", 22), rep("Zhang", 22), rep("PhenoAge", 22), rep("GrimAge", 22), rep("DunedinPACE", 22)), levels = c("Horvath", "Hannum", "Zhang", "PhenoAge", "GrimAge", "DunedinPACE"))
head(plot.data)

##       estimate  std.error     t.stat         p.val   ci.lower    ci.upper  PC
## PC1  0.6132778 0.08875804   6.909547  5.623812e-12  0.4393120  0.78724355 PC1
## PC2 -7.2567293 0.13546926 -53.567351  0.000000e+00 -7.5222490 -6.99120955 PC2
## PC3 -0.3002849 0.15930783  -1.884935  5.951021e-02 -0.6125283  0.01195843 PC3
## PC4  0.8487059 0.17677018   4.801183  1.634319e-06  0.5022364  1.19517544 PC4
## PC5 -2.0269365 0.20316377  -9.976860  3.557852e-23 -2.4251375 -1.62873548 PC5
## PC6 -6.8969005 0.20559171 -33.546589 6.572806e-218 -7.2998602 -6.49394076 PC6
##                     PC.label    type  clock
## PC1         PC1: Neutrophils DNAmAge Hannum
## PC2 PC2: T cell naive/memory DNAmAge Hannum
## PC3        PC3: Treg/CD4Tmem DNAmAge Hannum
## PC4           PC4: Bmem/Mono DNAmAge Hannum
## PC5             PC5: Bnv/Eos DNAmAge Hannum
## PC6   PC6: Naive CD8 T cells DNAmAge Hannum

ggplot(plot.data, aes(x = PC.label, y = estimate)) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_point(aes(color = type), shape = 19, size = 1.5) +
  geom_errorbar(aes(ymin = ci.lower, ymax = ci.upper, color = type), width = 0.4) +
  theme_bw() +
  scale_color_manual(values = c("#005580", "#993300"), drop = F) +
  scale_y_continuous(limits = c(-7.6, 4.7), breaks = seq(-10, 10, 2)) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        plot.margin = margin(t = 5.5, r = 5.5, b = 5.5, l = 22, unit = "pt")) +
  facet_wrap(facets = vars(clock), scales = "free", nrow = 2, dir = "h") +
  theme(panel.spacing.x = unit(1.2, "cm")) +
  labs(color= "", x = "", y = "")

## Warning: Removed 11 rows containing missing values (`geom_point()`).

---

Variance explained

Plot variance explained for each clock and residual.

plot.data <- var.df
plot.data$variable <- factor(plot.data$variable, levels = unique(plot.data$variable))

plot.data$group = factor(c("Age", "Age", "Hannum", "Hannum", "Horvath", "Horvath", "Zhang", "Zhang", "PhenoAge", "PhenoAge", "GrimAge", "GrimAge", "DunedinPACE", "DunedinPACE"), levels = c(c("Age", "Hannum", "Horvath", "Zhang", "PhenoAge", "GrimAge", "DunedinPACE")))
plot.data$type = factor(c("Age", "AgeAccel", "DNAmAge", "AgeAccel", "DNAmAge", "AgeAccel", "DNAmAge", "AgeAccel", "DNAmAge", "AgeAccel", "DNAmAge", "AgeAccel", "DNAmAge", "AgeAccel"), levels = c("Age", "DNAmAge", "AgeAccel"))
plot.data$var.round <- paste0(round(plot.data$var.ex, 0), "%")

#Split the data into 1st-gen and 2nd-gen clocks.
plot.data$gen <- c(rep("Age", 2), rep("1st-Generation Clocks", 6), rep("2nd-Generation Clocks", 6))
plot.data$gen <- factor(plot.data$gen, levels = c("Age", "1st-Generation Clocks", "2nd-Generation Clocks"))

#Remove unused rows.
plot.data$var.ex[c(2, 13)] <- 0
plot.data$var.round[c(2, 13)] <- ""
plot.data

##                         variable var.ex       group     type var.round
## age                          Age   43.5         Age      Age       44%
## ageRes                  AgeAccel    0.0         Age AgeAccel          
## hannum                    Hannum   52.6      Hannum  DNAmAge       53%
## hannumRes        Hannum residual   21.2      Hannum AgeAccel       21%
## horvath                  Horvath   43.8     Horvath  DNAmAge       44%
## horvathRes      Horvath residual    5.4     Horvath AgeAccel        5%
## zhang                      Zhang   46.0       Zhang  DNAmAge       46%
## zhangRes          Zhang residual    4.0       Zhang AgeAccel        4%
## phenoage                PhenoAge   49.7    PhenoAge  DNAmAge       50%
## phenoageRes    PhenoAge residual   19.9    PhenoAge AgeAccel       20%
## grimage                  GrimAge   46.6     GrimAge  DNAmAge       47%
## grimageRes      GrimAge residual   28.3     GrimAge AgeAccel       28%
## dunedinpaceAge   DunedinPACE-age    0.0 DunedinPACE  DNAmAge          
## dunedinpace          DunedinPACE   26.3 DunedinPACE AgeAccel       26%
##                                  gen
## age                              Age
## ageRes                           Age
## hannum         1st-Generation Clocks
## hannumRes      1st-Generation Clocks
## horvath        1st-Generation Clocks
## horvathRes     1st-Generation Clocks
## zhang          1st-Generation Clocks
## zhangRes       1st-Generation Clocks
## phenoage       2nd-Generation Clocks
## phenoageRes    2nd-Generation Clocks
## grimage        2nd-Generation Clocks
## grimageRes     2nd-Generation Clocks
## dunedinpaceAge 2nd-Generation Clocks
## dunedinpace    2nd-Generation Clocks

# Other visualization order.
plot.data <- plot.data[-c(2, 13),]

ggplot(data = plot.data, aes(x = group, y = var.ex, fill = type, label = var.round)) + 
  geom_bar(stat = "identity", width = 0.8, position = position_dodge(width = 0.8)) + 
  geom_text(hjust = 0.5, vjust = 1.15, position = position_dodge(width = 0.8), mapping = aes(color = type), size = 3) +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1), axis.ticks.x = element_blank(), panel.grid.major.x = element_blank()) +
  scale_fill_manual(values = c("#222222", "#005580", "#993300")) +
  scale_color_manual(values = c("white", "white", "white")) +
  scale_y_continuous(breaks = seq(0, 60, 10), limits = c(0, 53)) +
  labs(x = "", y = "Variance explained (%)", fill = "") +
  guides(color = "none", fill = "none") +
  facet_grid(cols = vars(type), scales = "free", space = "free")

Plot variance explained per PC.

PrintTopPCs <- function(){
  
  data.frame(
    row.names = paste0("PC", 1:11),
    var = rbind(
      round(summary(lm(formula = var ~ PC1, data = ss.pc))$r.squared * 100, 3),
      round(summary(lm(formula = var ~ PC2, data = ss.pc))$r.squared * 100, 3),
      round(summary(lm(formula = var ~ PC3, data = ss.pc))$r.squared * 100, 3),
      round(summary(lm(formula = var ~ PC3, data = ss.pc))$r.squared * 100, 3),
      round(summary(lm(formula = var ~ PC5, data = ss.pc))$r.squared * 100, 3),
      round(summary(lm(formula = var ~ PC6, data = ss.pc))$r.squared * 100, 3),
      round(summary(lm(formula = var ~ PC7, data = ss.pc))$r.squared * 100, 3),
      round(summary(lm(formula = var ~ PC8, data = ss.pc))$r.squared * 100, 3),
      round(summary(lm(formula = var ~ PC9, data = ss.pc))$r.squared * 100, 3),
      round(summary(lm(formula = var ~ PC10, data = ss.pc))$r.squared * 100, 3),
      round(summary(lm(formula = var ~ PC11, data = ss.pc))$r.squared * 100, 3)
      )
  )
  
}

top.pcs <- data.frame(
  row.names = paste0("PC", 1:11),
  PC = pc.labs,
  age = rep(NA, 11),
  hannum = rep(NA, 11),
  horvath = rep(NA, 11),
  zhang = rep(NA, 11),
  phenoage = rep(NA, 11),
  grimage = rep(NA, 11),
  dunedinpaceDummy = rep(NA, 11),
  ageRes = rep(NA, 11),
  hannumRes = rep(NA, 11),
  horvathRes = rep(NA, 11),
  zhangRes = rep(NA, 11),
  phenoageRes = rep(NA, 11),
  grimageRes = rep(NA, 11),
  dunedinpace = rep(NA, 11)
)

top.pcs$PC <- factor(top.pcs$PC, levels = (unique(top.pcs$PC)))

#Calendar age
ss.pc$var <- ss.pc$age
top.pcs$age <- PrintTopPCs()$var

#DNAmAge
ss.pc$var <- ss.pc$hannum
top.pcs$hannum <- PrintTopPCs()$var

ss.pc$var <- ss.pc$horvath
top.pcs$horvath <- PrintTopPCs()$var

ss.pc$var <- ss.pc$zhang
top.pcs$zhang <- PrintTopPCs()$var

ss.pc$var <- ss.pc$phenoage
top.pcs$phenoage <- PrintTopPCs()$var

ss.pc$var <- ss.pc$grimage
top.pcs$grimage <- PrintTopPCs()$var



#AgeAccels
ss.pc$var <- ss.pc$hannumRes
top.pcs$hannumRes <- PrintTopPCs()$var

ss.pc$var <- ss.pc$horvathRes
top.pcs$horvathRes <- PrintTopPCs()$var

ss.pc$var <- ss.pc$zhangRes
top.pcs$zhangRes <- PrintTopPCs()$var

ss.pc$var <- ss.pc$phenoageRes
top.pcs$phenoageRes <- PrintTopPCs()$var

ss.pc$var <- ss.pc$grimageRes
top.pcs$grimageRes <- PrintTopPCs()$var

ss.pc$var <- ss.pc$dunedinpace
top.pcs$dunedinpace <- PrintTopPCs()$var

#Remove unused columns.
top.pcs <- top.pcs[,-c(8, 9)]
top.pcs

##                               PC    age hannum horvath  zhang phenoage grimage
## PC1             PC1: Neutrophils  0.001  0.559   0.025  0.016    1.362   1.076
## PC2     PC2: T cell naive/memory 26.104 33.622  29.576 28.203   30.225  27.630
## PC3            PC3: Treg/CD4Tmem  0.240  0.042   0.110  0.183    0.025   0.230
## PC4               PC4: Bmem/Mono  0.240  0.042   0.110  0.183    0.025   0.230
## PC5                 PC5: Bnv/Eos  1.061  1.166   0.722  1.121    1.638   2.212
## PC6       PC6: Naive CD8 T cells 10.739 13.186   9.335 11.547   12.783  12.061
## PC7   PC7: Mixed (CD4Tnv/CD8Tnv)  4.725  3.353   3.242  4.317    3.193   2.488
## PC8                PC8: Eos/Mono  0.047  0.049   0.092  0.059    0.129   0.286
## PC9  PC9: Mixed (CD4Tnv+CD8Tmem)  0.140  0.100   0.263  0.174    0.014   0.000
## PC10          PC10: Mixed (Bmem)  0.001  0.024   0.002  0.001    0.001   0.004
## PC11          PC11: Mixed (Baso)  0.081  0.221   0.179  0.080    0.038   0.045
##      hannumRes horvathRes zhangRes phenoageRes grimageRes dunedinpace
## PC1      6.155      0.173    0.276      10.830     14.280       7.842
## PC2      9.800      3.499    2.427       4.139      1.553       8.211
## PC3      0.828      0.177    0.081       0.763      0.001       0.234
## PC4      0.828      0.177    0.081       0.763      0.001       0.234
## PC5      0.106      0.164    0.063       0.819      3.428       0.405
## PC6      2.910      0.031    0.913       2.093      1.389       7.590
## PC7      0.733      0.695    0.094       0.529      3.747       0.012
## PC8      0.003      0.096    0.027       0.202      1.496       1.127
## PC9      0.021      0.251    0.068       0.456      1.680       0.119
## PC10     0.168      0.001    0.001       0.000      0.131       0.000
## PC11     0.463      0.238    0.001       0.044      0.056       0.018

library(reshape2)
library(ggplot2)
plot.data <- melt(top.pcs, id.vars = "PC")
plot.data$group <- factor(NA, levels = c("Age", "DNAmAge",  "AgeAccel"))
plot.data$group[plot.data$variable %in% c("age", "ageRes")] <- "Age"
plot.data$group[plot.data$variable %in% c("hannum", "horvath", "zhang", "phenoage", "grimage", "dunedinpaceDummy")] <- "DNAmAge"
plot.data$group[plot.data$variable %in% c("hannumRes", "horvathRes", "zhangRes", "phenoageRes", "grimageRes", "dunedinpace")] <- "AgeAccel"

#Make a new variable that is the same for DNAmAge and AgeAccel of the same clock for better plotting.
plot.data$clock <- as.character(plot.data$variable)
plot.data$clock[grep("Res", plot.data$clock)] <- do.call("rbind", strsplit(plot.data$clock[grep("Res", plot.data$clock)], split = "Res"))
plot.data$clock[grep("Dummy", plot.data$clock)] <- do.call("rbind", strsplit(plot.data$clock[grep("Dummy", plot.data$clock)], split = "Dummy"))
plot.data$clock <- factor(plot.data$clock, levels = unique(plot.data$clock))
levels(plot.data$clock) <- c("Age", "Hannum", "Horvath", "Zhang", "PhenoAge", "GrimAge", "DunedinPACE")
head(plot.data)

##                         PC variable  value group clock
## 1         PC1: Neutrophils      age  0.001   Age   Age
## 2 PC2: T cell naive/memory      age 26.104   Age   Age
## 3        PC3: Treg/CD4Tmem      age  0.240   Age   Age
## 4           PC4: Bmem/Mono      age  0.240   Age   Age
## 5             PC5: Bnv/Eos      age  1.061   Age   Age
## 6   PC6: Naive CD8 T cells      age 10.739   Age   Age

ggplot(plot.data, aes(x = clock, y = PC, fill = value)) + 
  geom_tile() + 
  facet_grid(cols = vars(group), scales = "free", space = "free") +
  geom_text(aes(label = format(round(value, 1), 1)), vjust = 0.5, hjust = 0.5, size = 3.2) +
  scale_fill_gradient(low = "white", high = "red", limit = c(0, max(plot.data$value)),
                      guide = guide_colorbar(label.position = "left", label.hjust = 1)) +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        axis.title = element_blank(),
        panel.grid = element_blank(),
        strip.text = element_blank(),
        legend.title = element_text(hjust = 1),
        legend.position = "left") +
  scale_y_discrete(position = "right") +
  labs(x = "Clock variable", y = "Cell count PC", fill = "Variance\nexplained\n(%)")

---

SessionInfo

sessionInfo()

## R version 4.3.1 (2023-06-16)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Rocky Linux 8.10 (Green Obsidian)
## 
## Matrix products: default
## BLAS/LAPACK: /usr/lib64/libopenblas-r0.3.15.so;  LAPACK version 3.9.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: Europe/Amsterdam
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] ggpubr_0.6.0   reshape2_1.4.4 ggplot2_3.4.3 
## 
## loaded via a namespace (and not attached):
##  [1] sass_0.4.7        utf8_1.2.3        generics_0.1.3    tidyr_1.3.0      
##  [5] rstatix_0.7.2     stringi_1.7.12    lattice_0.21-8    digest_0.6.33    
##  [9] magrittr_2.0.3    evaluate_0.21     grid_4.3.1        fastmap_1.1.1    
## [13] plyr_1.8.8        jsonlite_1.8.7    Matrix_1.6-1      backports_1.4.1  
## [17] mgcv_1.8-42       purrr_1.0.2       fansi_1.0.4       scales_1.2.1     
## [21] jquerylib_0.1.4   abind_1.4-5       cli_3.6.1         rlang_1.1.1      
## [25] munsell_0.5.0     splines_4.3.1     withr_2.5.0       cachem_1.0.8     
## [29] yaml_2.3.7        tools_4.3.1       ggsignif_0.6.4    dplyr_1.1.2      
## [33] colorspace_2.1-0  broom_1.0.5       vctrs_0.6.3       R6_2.5.1         
## [37] lifecycle_1.0.3   stringr_1.5.0     car_3.1-2         pkgconfig_2.0.3  
## [41] pillar_1.9.0      bslib_0.5.1       gtable_0.3.4      glue_1.6.2       
## [45] Rcpp_1.0.11       highr_0.10        xfun_0.40         tibble_3.2.1     
## [49] tidyselect_1.2.0  rstudioapi_0.15.0 knitr_1.43        farver_2.1.1     
## [53] htmltools_0.5.6   nlme_3.1-162      labeling_0.4.2    rmarkdown_2.24   
## [57] carData_3.0-5     compiler_4.3.1